##################################
###   REPLICATION MATERIALS    ###
###            FOR             ###
### IMPLEMENTING PRESIDENTIAL  ###
###       PARTICULARISM        ###
##################################

# setwd("")
load("full_panel.RData")
load("aggregate.RData")
load("aggregate_nodefense.RData")
library(lfe)
library(ggplot2)
library(interflex)
library(gridExtra)
library(stargazer)
library(sandwich)
library(extrafont)
library(logr)

log_open(file_name = "replication_log.log")

plot.int <- function(obj, x1, x2, x1_string, x2_string, varcov){
  data <- model.frame(obj)
  attach(data, warn.conflicts = F)
  beta.hat <- coef(obj)
  z0 <- seq(floor(min(x2)), ceiling(max(x2)), length.out = 1000)
  dy.dx <- beta.hat[x1_string] + beta.hat[paste(x1_string,":",x2_string, sep="")]*z0
  se.dy.dx <- sqrt(varcov[x1_string, x1_string] + ((z0-mean(x2))^2)*varcov[x2_string, x2_string]
                   + 2*(abs(z0-mean(x2)))*varcov[x1_string, x2_string])
  upr <- dy.dx + 1.96*se.dy.dx
  lwr <- dy.dx - 1.96*se.dy.dx
  return(ggplot(data=NULL, aes(x=z0, y=dy.dx)) +
           geom_line(aes(z0, dy.dx),size = 1) +
           geom_hline(yintercept=0, size = 1, lty=2) +
           geom_ribbon(aes(ymin=lwr, ymax=upr), alpha=0.3) + theme_bw())
  detach(data)
}
predict.felm <- function(object, newdata, se.fit = FALSE,
                         interval = "none",
                         level = 0.95, vcov){
  if(missing(newdata)){
    stop("predict.felm requires newdata and predicts for all group effects = 0.")
  }
  
  tt <- terms(object)
  Terms <- delete.response(tt)
  attr(Terms, "intercept") <- 0
  
  m.mat <- model.matrix(Terms, data = newdata)
  m.coef <- as.numeric(object$coef)
  fit <- as.vector(m.mat %*% object$coef)
  fit <- data.frame(fit = fit)
  
  se.fit_mat <- sqrt(diag(m.mat %*%  vcov %*% t(m.mat)))
  if(interval == "confidence"){
    t_val <- qt((1 - level) / 2 + level, df = object$df.residual)
    fit$lwr <- fit$fit - t_val * se.fit_mat
    fit$upr <- fit$fit + t_val * se.fit_mat
  } else if (interval == "prediction"){
    stop("interval = \"prediction\" not yet implemented")
  }
  if(se.fit){
    return(list(fit=fit, se.fit=se.fit_mat))
  } else {
    return(fit)
  }
}

#### Table one ####

# Calculate median, within-agency standard deviations 
# of agency-president distance and politicization  
# Mummolo and Peterson (2019)

cj_dist_sd <- sd(felm(cj_dist ~ 0 |congress + agency_name | 0 | 0, data=panel_full)$residuals)
politicization_sd <- sd(felm(politicization ~ 0 |congress + agency_name | 0 | 0, data=panel_full)$residuals)

# Scale variables by median, within-agency standard deviation
panel_full$cj_dist_scaled <- NA
panel_full$politicization_scaled <- NA

panel_full$cj_dist_scaled <- 
  (panel_full$cj_dist - 
  mean(panel_full$cj_dist))/cj_dist_sd

panel_full$politicization_scaled <- 
  (panel_full$politicization - 
     mean(panel_full$politicization))/politicization_sd

agg1 <- felm(ln_funding ~ shared | name + congress | 0 | 0, aggregate)
agg2 <- felm(ln_funding ~ shared + majority + log(totalPopRaceFile) +
               log(medianIncome) + Appropriations + WaysAndMeans +
               close_election| name + congress | 0 | 0, aggregate)
mod1 <- felm(ln_funding ~ shared*cj_dist_scaled + shared*politicization_scaled|
               agency_name + congress| 0 | 0, data=subset(panel_full, placebo==0))
mod2 <- felm(ln_funding ~ shared*cj_dist_scaled + shared*politicization_scaled + majority +
               log(totalPopRaceFile) + log(medianIncome) + Appropriations + WaysAndMeans +
               close_election + pres_margin*shared + agency_mc_dist|
               agency_name + congress| 0 | agency_name, data=subset(panel_full, placebo==0))

vcov_agg1 <- vcovCL(agg1, aggregate$name)
vcov_agg2 <- vcovCL(agg2, aggregate$name)
vcov_mod1 <- vcovCL(mod1, panel_full$agency_name[panel_full$placebo==0])
vcov_mod2 <- vcovCL(mod2, panel_full$agency_name[panel_full$placebo==0])

ses <- lapply(list(vcov_agg1, vcov_agg2, vcov_mod1, vcov_mod2), diag)
ses <- lapply(ses, sqrt)

b2_b4_1 <- sum(mod1$coefficients[c(2,4)])
b3_b5_1 <- sum(mod1$coefficients[c(3,5)])
b2_b4_1_se <- sqrt(vcov_mod1[2,2] + vcov_mod1[4,4] + 2*vcov_mod1[2,4])
b3_b5_1_se <- sqrt(vcov_mod1[3,3] + vcov_mod1[5,5] + 2*vcov_mod1[3,5])

b2_b4_2 <- sum(mod2$coefficients[c(2,12)])
b3_b5_2 <- sum(mod2$coefficients[c(3,13)])
b2_b4_2_se <- sqrt(vcov_mod2[2,2] + vcov_mod2[12,12] + 2*vcov_mod2[2,12])
b3_b5_2_se <- sqrt(vcov_mod2[3,3] + vcov_mod2[13,13] + 2*vcov_mod2[3,13])

b2_b4_1_p <- 1-pnorm(abs(b2_b4_1/b2_b4_1_se))
b2_b4_2_p <- 1-pnorm(abs(b2_b4_2/b2_b4_2_se))
b3_b5_1_p <- 1-pnorm(abs(b3_b5_1/b3_b5_1_se))
b3_b5_2_p <- 1-pnorm(abs(b3_b5_2/b3_b5_2_se))

stargazer(agg1, agg2, mod1, mod2,
          keep = c("shared", "cj_dist_scaled", "politicization_scaled"),
          keep.stat = c("n", "adj.rsq"),
          covariate.labels = 
            c("Presidential Co-Partisan", "Agency--President Distance",
              "Politicization Ratio", "Pres. Co-Partisan \\times Ag.--Pres. Dist.",
              "Pres. Co-Partisan \\times Politicization Ratio"),
          se = ses,
          add.lines = list(c("beta2+beta4", "","",round(b2_b4_1, 3), round(b2_b4_2, 3)),
                           c("", "", "", paste("(",round(b2_b4_1_se, 3),")", sep=""), 
                                          paste("(",round(b2_b4_2_se,3), ")", sep="")),
                           c("beta3+beta5", "","",round(b3_b5_1, 3), round(b3_b5_2, 3)),
                           c("", "", "", paste("(",round(b3_b5_1_se, 3),")", sep=""), 
                             paste("(",round(b3_b5_2_se,3), ")", sep="")),
                           c("Congress FEs", "Yes", "Yes", "Yes", "Yes"),
                           c("Legislator FEs", "Yes", "Yes", "", ""),
                           c("Agency-Legislator FEs", "", "", "Yes", "Yes"),
                           c("Time-Varying Covariates", "", "Yes", "", "Yes")))

#### Figure one ####

cairo_pdf("marginal_effects_ideology.pdf", width=5, height=5)
plot.int(mod2, shared, cj_dist_scaled, "shared", "cj_dist_scaled",
         mod2$clustervcv) + 
  xlab("Agency-President Distance\n(Scaled by within-agency standard deviation)") +
  scale_x_continuous(expand = c(.01,.01)) +
  ylab("Marginal Effect of Presidential Co-Partisanship") +
  theme(panel.grid.minor = element_blank()) +
  geom_rug(aes(y=5,x=subset(panel_full, placebo==0)$cj_dist_scaled), size=2, sides="b",
           position = "jitter") + ylim(c(-2,2)) + xlim(c(-2,2)) +
  theme(text = element_text(family="Times New Roman", size = 14))
dev.off()

cairo_pdf("marginal_effects_politicization.pdf", width=5, height=5)
plot.int(mod2, shared, politicization_scaled, "shared", "politicization_scaled",
         vcovCL(mod2, panel_full$agency_name[panel_full$placebo==0])) + 
  xlab("Agency Politicization\n(Scaled by within-agency standard deviation)") +
  scale_x_continuous(expand = c(.01,.01)) +
  ylab("Marginal Effect of Presidential Co-Partisanship") +
  theme(panel.grid.minor = element_blank()) +
  geom_rug(aes(y=1,x=subset(panel_full, placebo==0)$politicization_scaled), size=2, sides="b",
           position = "jitter")+ ylim(c(-5,5)) +
  theme(text = element_text(family="Times New Roman", size = 14))
dev.off()


#### Table two ####
out <- data.frame(predict.felm(mod2,
                               data.frame(shared = rep(c(0,1), each=100),
                                          cj_dist_scaled = rep(seq(-2,2, length.out = 100), 2),
                                          totalPopRaceFile = mean(panel_full$totalPopRaceFile),
                                          medianIncome = mean(panel_full$medianIncome), majority =1,
                                          politicization_scaled = 0,
                                          agency_mc_dist = mean(panel_full$agency_mc_dist),
                                          Appropriations =0, WaysAndMeans =0, close_election=0,
                                          pres_margin = mean(panel_full$pres_margin),
                                          agency_name = ""), se.fit = T,
                               vcov = vcovCL(mod2,subset(panel_full, placebo==0)$agency_name)))

out$shared <- rep(c(0,1), each=100)
out$cj_dist_scaled <-rep(seq(-2,2, length.out = 100), 2)

sub_eff <- data.frame(cj_dist = c(-1, 0, 1),
                      co_part = c(exp(out$fit[125]+median(getfe(mod2, ef="zm2")$effect)-1)/1e+06,
                                  exp(out$fit[150]+median(getfe(mod2, ef="zm2")$effect)-1)/1e+06,
                                  exp(out$fit[175]+median(getfe(mod2, ef="zm2")$effect)-1)/1e+06),
                      contra_part = c(exp(out$fit[25]+median(getfe(mod2, ef="zm2")$effect)-1)/1e+06,
                                      exp(out$fit[50]+median(getfe(mod2, ef="zm2")$effect)-1)/1e+06,
                                      exp(out$fit[75]+median(getfe(mod2, ef="zm2")$effect)-1)/1e+06))
sub_eff$diff <- sub_eff$co_part - sub_eff$contra_part

sub_eff

#### Table three ####
mod1 <- felm(ln_funding ~ shared*cj_dist_scaled + shared*politicization_scaled|
               agency_name + congress| 0 | 0, data=subset(panel_full, placebo==0))
mod2 <- felm(ln_funding ~ shared*cj_dist_scaled + shared*politicization_scaled + majority +
               log(totalPopRaceFile) + log(medianIncome) + Appropriations + WaysAndMeans +
               close_election + pres_margin*shared + agency_mc_dist|
               agency_name + congress| 0 | 0, data=subset(panel_full, placebo==0))
mod3 <- felm(ln_funding ~ shared*cj_dist_scaled + shared*politicization_scaled|
               agency_name + congress| 0 | 0, data=subset(panel_full, placebo==1))
mod4 <- felm(ln_funding ~ shared*cj_dist_scaled + shared*politicization_scaled + majority +
               log(totalPopRaceFile) + log(medianIncome) + Appropriations + WaysAndMeans +
               close_election + pres_margin*shared + agency_mc_dist|
               agency_name + congress| 0 | 0, data=subset(panel_full, placebo==1))

vcov_mod3 <- vcovCL(mod3, panel_full$agency_name[panel_full$placebo==1])
vcov_mod4 <- vcovCL(mod4, panel_full$agency_name[panel_full$placebo==1])

vcovs <- list(
  vcov_mod1, vcov_mod2,
  vcov_mod3, vcov_mod4
)

ses <- lapply(vcovs, diag)
ses <- lapply(ses, sqrt)

b2_b4_3 <- sum(mod3$coefficients[c(2,4)])
b3_b5_3 <- sum(mod3$coefficients[c(3,5)])
b2_b4_3_se <- sqrt(vcov_mod3[2,2] + vcov_mod3[4,4] + 2*vcov_mod3[2,4])
b3_b5_3_se <- sqrt(vcov_mod3[3,3] + vcov_mod3[5,5] + 2*vcov_mod3[3,5])

b2_b4_4 <- sum(mod4$coefficients[c(2,12)])
b3_b5_4 <- sum(mod4$coefficients[c(3,13)])
b2_b4_4_se <- sqrt(vcov_mod4[2,2] + vcov_mod4[12,12] + 2*vcov_mod4[2,12])
b3_b5_4_se <- sqrt(vcov_mod4[3,3] + vcov_mod4[13,13] + 2*vcov_mod4[3,13])

b2_b4_3_p <- 1-pnorm(abs(b2_b4_3/b2_b4_3_se))
b2_b4_4_p <- 1-pnorm(abs(b2_b4_4/b2_b4_4_se))
b3_b5_3_p <- 1-pnorm(abs(b3_b5_3/b3_b5_3_se))
b3_b5_4_p <- 1-pnorm(abs(b3_b5_4/b3_b5_4_se))

stargazer(mod1, mod2, mod3, mod4,
          keep = c("shared", "cj_dist_scaled", "politicization_scaled"),
          keep.stat = c("n", "adj.rsq"),
          covariate.labels = 
            c("Presidential Co-Partisan", "Agency--President Distance",
              "Politicization Ratio", "Pres. Co-Partisan \\times Ag.--Pres. Dist.",
              "Pres. Co-Partisan \\times Politicization Ratio"),
          se = ses,
          add.lines = list(c("beta2+beta4",round(b2_b4_1, 3), round(b2_b4_2, 3),
                             round(b2_b4_3, 3), round(b2_b4_4, 3)),
                           c("", paste("(",round(b2_b4_1_se, 3),")", sep=""), 
                             paste("(",round(b2_b4_2_se,3), ")", sep=""),
                             paste("(",round(b2_b4_3_se, 3),")", sep=""),
                             paste("(",round(b2_b4_4_se, 3),")", sep="")),
                           c("beta3+beta5",round(b3_b5_1, 3), round(b3_b5_2, 3),
                             round(b3_b5_3, 3), round(b3_b5_2, 4)),
                           c("", paste("(",round(b3_b5_1_se, 3),")", sep=""), 
                             paste("(",round(b3_b5_2_se,3), ")", sep=""),
                             paste("(",round(b3_b5_3_se, 3),")", sep=""), 
                             paste("(",round(b3_b5_4_se,3), ")", sep="")),
                           c("Congress FEs", rep("Yes", 4)),
                           c("Agency-Legislator FEs", rep("Yes", 4)),
                           c("Time-Varying Covariates", "", "Yes", "", "Yes")))


#### Figure two ####

cairo_pdf("marginal_effects_ideology_placebo.pdf", width=5, height=5)
plot.int(mod4, shared, cj_dist_scaled, "shared", "cj_dist_scaled",
         vcovCL(mod4, panel_full$agency_name[panel_full$placebo==0])) + 
  xlab("Agency-President Distance\n(Scaled by within-agency standard deviation)") +
  scale_x_continuous(breaks=seq(0,1,length.out = 9),labels = format(seq(-2,2,.5), nsmall = 1),
                     expand = c(.01,.01)) +
  ylab("Marginal Effect of Presidential Co-Partisanship") +
  theme(panel.grid.minor = element_blank()) +
  geom_rug(aes(y=1,x=subset(panel_full, placebo==1)$cj_dist_scaled), size=2, sides="b",
           position = "jitter")+ ylim(c(-2,2)) + xlim(c(-2,2)) +
  theme(text = element_text(family="Times New Roman", size = 14))
dev.off()

cairo_pdf("marginal_effects_politicization_placebo.pdf", width=5, height=5)
plot.int(mod4, shared, politicization_scaled, "shared", "politicization_scaled",
         vcovCL(mod4, panel_full$agency_name[panel_full$placebo==0])) + 
  xlab("Agency Politicization\n(Scaled by within-agency standard deviation)") +
  scale_x_continuous(expand = c(.01,.01)) +
  ylab("Marginal Effect of Presidential Co-Partisanship") +
  theme(panel.grid.minor = element_blank()) +
  geom_rug(aes(y=1,x=subset(panel_full, placebo==1)$politicization_scaled), size=2, sides="b",
           position = "jitter")+ ylim(c(-5,5)) +
  theme(text = element_text(family="Times New Roman", size = 14))
dev.off()


##################
#### APPENDIX ####
##################

#### Interflex ####

panel_full$shared_polit <- panel_full$shared*panel_full$politicization_scaled
panel_full$shared_cj <- panel_full$shared*panel_full$cj_dist_scaled
panel_full$shared_marg <- panel_full$shared*panel_full$pres_marg
panel_full$ln_pop <- log(panel_full$totalPopRaceFile)
panel_full$ln_income <- log(panel_full$medianIncome)

out <- inter.binning(data=subset(panel_full, placebo==0),
                     Y="ln_funding", D="shared", X="cj_dist_scaled",
                     Z=c("majority", "politicization_scaled",
                         "agency_mc_dist", "Appropriations",
                         "WaysAndMeans", "close_election",
                         "pres_margin", "ln_pop", "ln_income",
                         "shared_polit", "shared_marg"),
                     FE=c("agency_name", "congress"),
                     vartype="cluster", cl="agency_name",
                     nbins = 5)
cairo_pdf("interflex.pdf", width=10, height=5)
plot(out, xlab = "Agency-President Distance",
     ylab = "Marginal Effect of Presidential Co-Partisanship")
dev.off()

#### FE Adjustments ####

# Ideology #
FEs <- felm(cj_dist ~ 0 |congress + agency_name | 0 | 0, data=panel_full)
FEadj <- data.frame(resids = FEs$residuals, unadjust = scale(panel_full$cj_dist, center=T, scale=F))
names(FEadj)[1] <- "resids"

pdf("FEadj_ideology.pdf", width=5, height=5)
ggplot(FEadj) + geom_density(aes(x=resids), fill="red", alpha=.5) +
  geom_density(aes(x=unadjust), fill="grey", alpha=.5) + theme_bw() + 
  xlab("Agency-President Distance \n (Centered on Zero)") + ylab("Density") +
  annotate(.15, 5, label="After FEs \n SD = 0.102", color="red", geom="text") + 
  annotate(-.35, 1.75, label="Original Distribution\nSD = 0.263", geom="text") +
  theme(panel.grid = element_blank())
dev.off()

panel_full$agency <- as.character(panel_full$agency)
within <- data.frame(within=tapply(panel_full$cj_dist, panel_full$agency, max)-tapply(panel_full$cj_dist, panel_full$agency, min))
pdf("within_ideology.pdf", width=5, height=5)
ggplot(within) + geom_histogram(aes(x=within), fill="grey", bins=10, col="black") +
  theme_bw() + geom_vline(xintercept = median(within$within), lty=2) + ylim(c(0,5)) +
  geom_vline(xintercept = quantile(within$within, .95), lty=2) +
  geom_segment(aes(x=.51, xend=median(within), y=4,yend=4),
               arrow=arrow(angle=30, ends = "last", length=unit(.1, "in"))) + 
  annotate(geom="text", .55, 4.15, label="Median") + 
  geom_segment(aes(x=.77, xend=quantile(within, .95), y=2.5,yend=2.5),
               arrow=arrow(angle=30, ends = "last", length=unit(.1, "in"))) + 
  annotate(geom="text", .815, 2.8, label="95th\nPercentile") + 
  xlab("Within-Agency Range in Agency-President Distance") + ylab("Count")
dev.off()


FEs <- felm(politicization ~ 0 |congress + agency_name | 0 | 0, data=panel_full)
FEadj <- data.frame(resids = FEs$residuals, unadjust = scale(panel_full$politicization, center=T, scale=F))
names(FEadj)[1] <- "resids"

pdf("FEadj_polticization.pdf", width=5, height=5)
ggplot(FEadj) + geom_density(aes(x=resids), fill="red", alpha=.5) +
  geom_density(aes(x=unadjust), fill="grey", alpha=.5) + theme_bw() + 
  xlab("Agency Politicization \n (Centered on Zero)") + ylab("Density") +
  annotate(.25, 5, label="After FEs \n SD = 0.088", color="red", geom="text") + 
  annotate(-.45, 2.5, label="Original Distribution\nSD = 0.494", geom="text") +
  theme(panel.grid = element_blank()) + xlim(c(-1,3))
dev.off()

panel_full$agency <- as.character(panel_full$agency)
within <- data.frame(within=tapply(panel_full$politicization, panel_full$agency, max)-tapply(panel_full$politicization, panel_full$agency, min))
pdf("within_politicization.pdf", width=5, height=5)
ggplot(within) + geom_histogram(aes(x=within), fill="grey", bins=10, col="black") +
  theme_bw() + geom_vline(xintercept = median(within$within), lty=2) +
  geom_vline(xintercept = quantile(within$within, .95), lty=2) +
  geom_segment(aes(x=.51, xend=median(within), y=4,yend=4),
               arrow=arrow(angle=30, ends = "last", length=unit(.1, "in"))) + 
  annotate(geom="text", .55, 4.15, label="Median") + 
  geom_segment(aes(x=.77, xend=quantile(within, .95), y=2.5,yend=2.5),
               arrow=arrow(angle=30, ends = "last", length=unit(.1, "in"))) + 
  annotate(geom="text", .815, 2.8, label="95th\nPercentile") + 
  xlab("Within-Agency Range in Agency Politicization") + ylab("Count")
dev.off()

#### Alternative Clusters ####
mod1 <- felm(ln_funding ~ shared*cj_dist_scaled + shared*politicization_scaled|
               agency_name + congress| 0 | 0, data=subset(panel_full, placebo==0))
mod2 <- felm(ln_funding ~ shared*cj_dist_scaled + shared*politicization_scaled + majority +
               log(totalPopRaceFile) + log(medianIncome) + Appropriations + WaysAndMeans +
               close_election + pres_margin*shared + agency_mc_dist|
               agency_name + congress| 0 | agency_name, data=subset(panel_full, placebo==0))

vcov_mod1 <- vcovCL(mod1, panel_full$agency[panel_full$placebo==0])
vcov_mod2 <- vcovCL(mod2, panel_full$agency[panel_full$placebo==0])
vcov_mod3 <- vcovCL(mod1, panel_full$name[panel_full$placebo==0])
vcov_mod4 <- vcovCL(mod2, panel_full$name[panel_full$placebo==0])

ses <- lapply(list(vcov_mod1, vcov_mod2, vcov_mod3, vcov_mod4), diag)
ses <- lapply(ses, sqrt)

b2_b4_1 <- sum(mod1$coefficients[c(2,4)])
b3_b5_1 <- sum(mod1$coefficients[c(3,5)])
b2_b4_1_se <- sqrt(vcov_mod1[2,2] + vcov_mod1[4,4] + 2*vcov_mod1[2,4])
b3_b5_1_se <- sqrt(vcov_mod1[3,3] + vcov_mod1[3,5] + 2*vcov_mod1[3,5])

b2_b4_2 <- sum(mod2$coefficients[c(2,12)])
b3_b5_2 <- sum(mod2$coefficients[c(3,13)])
b2_b4_2_se <- sqrt(vcov_mod2[2,2] + vcov_mod2[12,12] + 2*vcov_mod2[2,12])
b3_b5_2_se <- sqrt(vcov_mod2[3,3] + vcov_mod2[13,13] + 2*vcov_mod2[3,13])


b2_b4_3 <- sum(mod1$coefficients[c(2,4)])
b3_b5_3 <- sum(mod1$coefficients[c(3,5)])
b2_b4_3_se <- sqrt(vcov_mod3[2,2] + vcov_mod3[4,4] + 2*vcov_mod3[2,4])
b3_b5_3_se <- sqrt(vcov_mod3[3,3] + vcov_mod3[5,5] + 2*vcov_mod3[3,5])

b2_b4_4 <- sum(mod2$coefficients[c(2,12)])
b3_b5_4 <- sum(mod2$coefficients[c(3,13)])
b2_b4_4_se <- sqrt(vcov_mod4[2,2] + vcov_mod4[12,12] + 2*vcov_mod4[2,12])
b3_b5_4_se <- sqrt(vcov_mod4[3,3] + vcov_mod4[13,13] + 2*vcov_mod4[3,13])

b2_b4_1_p <- 1-pnorm(abs(b2_b4_1/b2_b4_1_se))
b2_b4_2_p <- 1-pnorm(abs(b2_b4_2/b2_b4_2_se))
b3_b5_1_p <- 1-pnorm(abs(b3_b5_1/b3_b5_1_se))
b3_b5_2_p <- 1-pnorm(abs(b3_b5_2/b3_b5_2_se))
b2_b4_3_p <- 1-pnorm(abs(b2_b4_3/b2_b4_3_se))
b2_b4_4_p <- 1-pnorm(abs(b2_b4_4/b2_b4_4_se))
b3_b5_3_p <- 1-pnorm(abs(b3_b5_3/b3_b5_3_se))
b3_b5_4_p <- 1-pnorm(abs(b3_b5_4/b3_b5_4_se))

stargazer(mod1, mod2, mod1, mod2,
          keep = c("shared", "cj_dist_scaled", "politicization_scaled"),
          keep.stat = c("n", "adj.rsq"),
          covariate.labels = 
            c("Presidential Co-Partisan", "Agency--President Distance",
              "Politicization Ratio", "Pres. Co-Partisan \\times Ag.--Pres. Dist.",
              "Pres. Co-Partisan \\times Politicization Ratio"),
          se = ses,
          add.lines = list(c("beta2+beta4",round(b2_b4_1, 3), round(b2_b4_2, 3),
                             round(b2_b4_3, 3), round(b2_b4_4, 3)),
                           c("", paste("(",round(b2_b4_1_se, 3),")", sep=""), 
                             paste("(",round(b2_b4_2_se,3), ")", sep=""),
                             paste("(",round(b2_b4_3_se, 3),")", sep=""), 
                             paste("(",round(b2_b4_4_se,3), ")", sep="")),
                           c("beta3+beta5", round(b3_b5_1, 3), round(b3_b5_2, 3),
                             round(b3_b5_3, 3), round(b3_b5_4, 3)),
                           c("", paste("(",round(b3_b5_1_se, 3),")", sep=""), 
                             paste("(",round(b3_b5_2_se,3), ")", sep=""),
                             paste("(",round(b3_b5_3_se, 3),")", sep=""), 
                             paste("(",round(b3_b5_4_se,3), ")", sep="")),
                           c("Congress FEs", rep("Yes", 4)),
                           c("Agency-Legislator FEs", rep("Yes", 4)),
                           c("Time-Varying Covariates", "", "Yes", "", "Yes")))


#### High-Variance Programs ####
mod1 <- felm(ln_highvar_funding ~ shared*cj_dist_scaled + shared*politicization_scaled|
               agency_name + congress| 0 | 0, data=subset(panel_full, placebo==0))
mod2 <- felm(ln_highvar_funding ~ shared*cj_dist_scaled + shared*politicization_scaled + majority +
               log(totalPopRaceFile) + log(medianIncome) + Appropriations + WaysAndMeans +
               close_election + pres_margin*shared + agency_mc_dist|
               agency_name + congress| 0 | agency_name, data=subset(panel_full, placebo==0))

vcov_mod1 <- vcovCL(mod1, panel_full$agency_name[panel_full$placebo==0])
vcov_mod2 <- vcovCL(mod2, panel_full$agency_name[panel_full$placebo==0])

vcovs <- list(
  vcov_mod1, vcov_mod2
)

ses <- lapply(vcovs, diag)
ses <- lapply(ses, sqrt)

b2_b4_1 <- sum(mod1$coefficients[c(2,4)])
b3_b5_1 <- sum(mod1$coefficients[c(3,5)])
b2_b4_1_se <- sqrt(vcov_mod1[2,2] + vcov_mod1[4,4] + 2*vcov_mod1[2,4])
b3_b5_1_se <- sqrt(vcov_mod1[3,3] + vcov_mod1[5,5] + 2*vcov_mod1[3,5])

b2_b4_2 <- sum(mod2$coefficients[c(2,12)])
b3_b5_2 <- sum(mod2$coefficients[c(3,13)])
b2_b4_2_se <- sqrt(vcov_mod2[2,2] + vcov_mod2[12,12] + 2*vcov_mod2[2,12])
b3_b5_2_se <- sqrt(vcov_mod2[3,3] + vcov_mod2[13,13] + 2*vcov_mod2[3,13])

b2_b4_1_p <- 1-pnorm(abs(b2_b4_1/b2_b4_1_se))
b2_b4_2_p <- 1-pnorm(abs(b2_b4_2/b2_b4_2_se))
b3_b5_1_p <- 1-pnorm(abs(b3_b5_1/b3_b5_1_se))
b3_b5_2_p <- 1-pnorm(abs(b3_b5_2/b3_b5_2_se))

stargazer(mod1, mod2, 
          keep = c("shared", "cj_dist_scaled", "politicization_scaled"),
          keep.stat = c("n", "adj.rsq"),
          covariate.labels = 
            c("Presidential Co-Partisan", "Agency--President Distance",
              "Politicization Ratio", "Pres. Co-Partisan \\times Ag.--Pres. Dist.",
              "Pres. Co-Partisan \\times Politicization Ratio"),
          se = ses,
          add.lines = list(c("beta2+beta4",round(b2_b4_1, 3), round(b2_b4_2, 3)),
                           c("", paste("(",round(b2_b4_1_se, 3),")", sep=""), 
                             paste("(",round(b2_b4_2_se,3), ")", sep="")),
                           c("beta3+beta5", 
                             round(b3_b5_3, 3), round(b3_b5_4, 3)),
                           c("", paste("(",round(b3_b5_1_se, 3),")", sep=""), 
                             paste("(",round(b3_b5_2_se,3), ")", sep="")),
                           c("Congress FEs", rep("Yes", 2)),
                           c("Agency-Legislator FEs", rep("Yes", 2)),
                           c("Time-Varying Covariates", "", "Yes")))

#### Excluding Defense Agencies ####
panel_nodefense <- subset(panel_full,
                          !(agency %in% c("Department of Defense", "Army", "Navy", "Air Force")))


mod1 <- felm(ln_funding ~ shared*cj_dist_scaled + shared*politicization_scaled|
               agency_name + congress| 0 | 0, data=subset(panel_nodefense, placebo==0))
mod2 <- felm(ln_funding ~ shared*cj_dist_scaled + shared*politicization_scaled + majority +
               log(totalPopRaceFile) + log(medianIncome) + Appropriations + WaysAndMeans +
               close_election + pres_margin*shared + agency_mc_dist|
               agency_name + congress| 0 | 0, data=subset(panel_nodefense, placebo==0))
mod3 <- felm(ln_funding ~ shared*cj_dist_scaled + shared*politicization_scaled|
               agency_name + congress| 0 | 0, data=subset(panel_nodefense, placebo==1))
mod4 <- felm(ln_funding ~ shared*cj_dist_scaled + shared*politicization_scaled + majority +
               log(totalPopRaceFile) + log(medianIncome) + Appropriations + WaysAndMeans +
               close_election + pres_margin*shared + agency_mc_dist|
               agency_name + congress| 0 | 0, data=subset(panel_nodefense, placebo==1))

vcov_mod1 <- vcovCL(mod1, panel_nodefense$agency_name[panel_nodefense$placebo==0])
vcov_mod2 <- vcovCL(mod2, panel_nodefense$agency_name[panel_nodefense$placebo==0])
vcov_mod3 <- vcovCL(mod3, panel_nodefense$agency_name[panel_nodefense$placebo==1])
vcov_mod4 <- vcovCL(mod4, panel_nodefense$agency_name[panel_nodefense$placebo==1])

vcovs <- list(
  vcov_mod1, vcov_mod2,
  vcov_mod3, vcov_mod4
)

ses <- lapply(vcovs, diag)
ses <- lapply(ses, sqrt)

b2_b4_1 <- sum(mod1$coefficients[c(2,4)])
b3_b5_1 <- sum(mod1$coefficients[c(3,5)])
b2_b4_1_se <- sqrt(vcov_mod1[2,2] + vcov_mod1[4,4] + 2*vcov_mod1[2,4])
b3_b5_1_se <- sqrt(vcov_mod1[3,3] + vcov_mod1[3,5] + 2*vcov_mod1[3,5])

b2_b4_2 <- sum(mod2$coefficients[c(2,12)])
b3_b5_2 <- sum(mod2$coefficients[c(3,13)])
b2_b4_2_se <- sqrt(vcov_mod2[2,2] + vcov_mod2[12,12] + 2*vcov_mod2[2,12])
b3_b5_2_se <- sqrt(vcov_mod2[3,3] + vcov_mod2[13,13] + 2*vcov_mod2[3,13])

b2_b4_1_p <- 1-pnorm(abs(b2_b4_1/b2_b4_1_se))
b2_b4_2_p <- 1-pnorm(abs(b2_b4_2/b2_b4_2_se))
b3_b5_1_p <- 1-pnorm(abs(b3_b5_1/b3_b5_1_se))
b3_b5_2_p <- 1-pnorm(abs(b3_b5_2/b3_b5_2_se))

b2_b4_3 <- sum(mod3$coefficients[c(2,4)])
b3_b5_3 <- sum(mod3$coefficients[c(3,5)])
b2_b4_3_se <- sqrt(vcov_mod3[2,2] + vcov_mod3[4,4] + 2*vcov_mod3[2,4])
b3_b5_3_se <- sqrt(vcov_mod3[3,3] + vcov_mod3[5,5] + 2*vcov_mod3[3,5])

b2_b4_4 <- sum(mod4$coefficients[c(2,12)])
b3_b5_4 <- sum(mod4$coefficients[c(3,13)])
b2_b4_4_se <- sqrt(vcov_mod4[2,2] + vcov_mod4[12,12] + 2*vcov_mod4[2,12])
b3_b5_4_se <- sqrt(vcov_mod4[3,3] + vcov_mod4[13,13] + 2*vcov_mod4[3,13])

b2_b4_3_p <- 1-pnorm(abs(b2_b4_3/b2_b4_3_se))
b2_b4_4_p <- 1-pnorm(abs(b2_b4_4/b2_b4_4_se))
b3_b5_3_p <- 1-pnorm(abs(b3_b5_3/b3_b5_3_se))
b3_b5_4_p <- 1-pnorm(abs(b3_b5_4/b3_b5_4_se))

stargazer(mod1, mod2, mod3, mod4,
          keep = c("shared", "cj_dist_scaled", "politicization_scaled"),
          keep.stat = c("n", "adj.rsq"),
          covariate.labels = 
            c("Presidential Co-Partisan", "Agency--President Distance",
              "Politicization Ratio", "Pres. Co-Partisan \\times Ag.--Pres. Dist.",
              "Pres. Co-Partisan \\times Politicization Ratio"),
          se = ses,
          add.lines = list(c("beta2+beta4",round(b2_b4_1, 3), round(b2_b4_2, 3),
                             round(b2_b4_3, 3), round(b2_b4_4, 3)),
                           c("", paste("(",round(b2_b4_1_se, 3),")", sep=""), 
                             paste("(",round(b2_b4_2_se,3), ")", sep=""),
                             paste("(",round(b2_b4_3_se, 3),")", sep=""),
                             paste("(",round(b2_b4_4_se, 3),")", sep="")),
                           c("beta3+beta5",round(b3_b5_1, 3), round(b3_b5_2, 3),
                             round(b3_b5_3, 3), round(b3_b5_2, 4)),
                           c("", paste("(",round(b3_b5_1_se, 3),")", sep=""), 
                             paste("(",round(b3_b5_2_se,3), ")", sep=""),
                             paste("(",round(b3_b5_3_se, 3),")", sep=""), 
                             paste("(",round(b3_b5_4_se,3), ")", sep="")),
                           c("Congress FEs", rep("Yes", 4)),
                           c("Agency-Legislator FEs", rep("Yes", 4)),
                           c("Time-Varying Covariates", "", "Yes", "", "Yes")))

#### By Structure ####
panel_cabinet <- panel_full[grepl("Department", panel_full$agency),]
panel_ind <- panel_full[!grepl("Department", panel_full$agency),]

mod1 <- felm(ln_funding ~ shared*cj_dist_scaled + shared*politicization_scaled|
               agency_name + congress| 0 | 0, data=subset(panel_cabinet, placebo==0))
mod2 <- felm(ln_funding ~ shared*cj_dist_scaled + shared*politicization_scaled + majority +
               log(totalPopRaceFile) + log(medianIncome) + Appropriations + WaysAndMeans +
               close_election + pres_margin*shared + agency_mc_dist|
               agency_name + congress| 0 | agency_name, data=subset(panel_cabinet, placebo==0))
mod3 <- felm(ln_funding ~ shared*cj_dist_scaled + shared*politicization_scaled|
               agency_name + congress| 0 | 0, data=subset(panel_ind, placebo==0))
mod4 <- felm(ln_funding ~ shared*cj_dist_scaled + shared*politicization_scaled + majority +
               log(totalPopRaceFile) + log(medianIncome) + Appropriations + WaysAndMeans +
               close_election + pres_margin*shared + agency_mc_dist|
               agency_name + congress| 0 | agency_name, data=subset(panel_ind, placebo==0))

vcov_mod1 <- vcovCL(mod1, panel_cabinet$agency[panel_cabinet$placebo==0])
vcov_mod2 <- vcovCL(mod2, panel_cabinet$agency[panel_cabinet$placebo==0])
vcov_mod3 <- vcovCL(mod3, panel_ind$name[panel_ind$placebo==0])
vcov_mod4 <- vcovCL(mod4, panel_ind$name[panel_ind$placebo==0])

ses <- lapply(list(vcov_mod1, vcov_mod2, vcov_mod3, vcov_mod4), diag)
ses <- lapply(ses, sqrt)

b2_b4_1 <- sum(mod1$coefficients[c(2,4)])
b3_b5_1 <- sum(mod1$coefficients[c(3,5)])
b2_b4_1_se <- sqrt(vcov_mod1[2,2] + vcov_mod1[4,4] + 2*vcov_mod1[2,4])
b3_b5_1_se <- sqrt(vcov_mod1[3,3] + vcov_mod1[5,5] + 2*vcov_mod1[3,5])

b2_b4_2 <- sum(mod2$coefficients[c(2,12)])
b3_b5_2 <- sum(mod2$coefficients[c(3,13)])
b2_b4_2_se <- sqrt(vcov_mod2[2,2] + vcov_mod2[12,12] + 2*vcov_mod2[2,12])
b3_b5_2_se <- sqrt(vcov_mod2[3,3] + vcov_mod2[13,13] + 2*vcov_mod2[3,13])


b2_b4_3 <- sum(mod3$coefficients[c(2,4)])
b3_b5_3 <- sum(mod3$coefficients[c(3,5)])
b2_b4_3_se <- sqrt(vcov_mod3[2,2] + vcov_mod3[4,4] + 2*vcov_mod3[2,4])
b3_b5_3_se <- sqrt(vcov_mod3[3,3] + vcov_mod3[5,5] + 2*vcov_mod3[3,5])

b2_b4_4 <- sum(mod4$coefficients[c(2,12)])
b3_b5_4 <- sum(mod4$coefficients[c(3,13)])
b2_b4_4_se <- sqrt(vcov_mod4[2,2] + vcov_mod4[12,12] + 2*vcov_mod4[2,12])
b3_b5_4_se <- sqrt(vcov_mod4[3,3] + vcov_mod4[13,13] + 2*vcov_mod4[3,13])

b2_b4_1_p <- 1-pnorm(abs(b2_b4_1/b2_b4_1_se))
b2_b4_2_p <- 1-pnorm(abs(b2_b4_2/b2_b4_2_se))
b3_b5_1_p <- 1-pnorm(abs(b3_b5_1/b3_b5_1_se))
b3_b5_2_p <- 1-pnorm(abs(b3_b5_2/b3_b5_2_se))
b2_b4_3_p <- 1-pnorm(abs(b2_b4_3/b2_b4_3_se))
b2_b4_4_p <- 1-pnorm(abs(b2_b4_4/b2_b4_4_se))
b3_b5_3_p <- 1-pnorm(abs(b3_b5_3/b3_b5_3_se))
b3_b5_4_p <- 1-pnorm(abs(b3_b5_4/b3_b5_4_se))

stargazer(mod1, mod2, mod3, mod4, 
          keep = c("shared", "cj_dist_scaled", "politicization_scaled"),
          keep.stat = c("n", "adj.rsq"),
          covariate.labels = 
            c("Presidential Co-Partisan", "Agency--President Distance",
              "Politicization Ratio", "Pres. Co-Partisan \\times Ag.--Pres. Dist.",
              "Pres. Co-Partisan \\times Politicization Ratio"),
          se = ses,
          add.lines = list(c("beta2+beta4",round(b2_b4_1, 3), round(b2_b4_2, 3),
                             round(b2_b4_3, 3), round(b2_b4_4, 3)),
                           c("", paste("(",round(b2_b4_1_se, 3),")", sep=""), 
                             paste("(",round(b2_b4_2_se,3), ")", sep=""),
                             paste("(",round(b2_b4_3_se, 3),")", sep=""), 
                             paste("(",round(b2_b4_4_se,3), ")", sep="")),
                           c("beta3+beta5", round(b3_b5_1, 3), round(b3_b5_2, 3),
                             round(b3_b5_3, 3), round(b3_b5_4, 3)),
                           c("", paste("(",round(b3_b5_1_se, 3),")", sep=""), 
                             paste("(",round(b3_b5_2_se,3), ")", sep=""),
                             paste("(",round(b3_b5_3_se, 3),")", sep=""), 
                             paste("(",round(b3_b5_4_se,3), ")", sep="")),
                           c("Congress FEs", rep("Yes", 4)),
                           c("Agency-Legislator FEs", rep("Yes", 4)),
                           c("Time-Varying Covariates", "", "Yes", "", "Yes")))

log_close()

